Calculation of nonzero-temperature Casimir forces in the time domain 
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We show how to compute Casimir forces at nonzero temperatures with time-domain electromag- 
netic simulations, for example using a finite- difference time-domain (FDTD) method. Compared 
to our previous zero-temperature time-domain method, only a small modification is required, but 
we explain that some care is required to properly capture the zero-frequency contribution. We 
validate the method against analytical and numerical frequency-domain calculations, and show a 
surprising high-temperature disappearance of a non-monotonic behavior previously demonstrated 
in a piston-like geometry. 



In this paper, we show how to compute nonzero- 
temperature (T > 0) corrections to Casimir forces via 
time-domain calculations, generalizing a computational 
approach based on the finite-difference time-domain 
(FDTD) method that we previously demonstrated for 
T = New computational methods for Casimir in- 

teractions have become important in order to model non- 
planar micromechanical systems where unusual Casimir 
effects have been predicted [ll-[l2j, and there has been 
increasing interest in T > corrections [l^4l9j|. espe- 
cially in recently identified systems where these effects 
are non- negligible [l6j|. Although T > effects are easy 
to incorporate in the imaginary frequency domain, where 
they merely turn an integral into a sum over Matsubara 
frequencies [20[, they turn out to be nontrivial to han- 
dle in the time-domain because of the singularity of the 
zero-frequency contribution, and we show that a naive 
approach leads to incorrect results. We validate our ap- 
proach both with a one-dimensional system where analyt- 
ical solutions are available, and also in a two-dimensional 
(2D) piston-like geometry [3, H, Hl|, [22| where we com- 
pare to a frequency-domain numerical method. In the 
2D piston geometry, we observe an interesting effect in 
which a non-monotonic phenomenon previously identi- 
fied at T = disappears for a sufficiently large T. 

The Casimir force is a combination of fluctuations at 
all frequencies cj, and the T = force can be expressed as 
an integral F(0) = J °° f(t;)d£ over Wick-rotated imag- 
inary frequencies uj = z£ [20j. At a nonzero T, this 
integral is replaced by a sum over "Matsubara frequen- 
cies" £ n = uttujt for integers n, where ujt = 2k^T /h and 
ks is Boltzmann's constant [201 ] : 



F(T) = ttujt 
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The transformation from the T = integral to a sum- 
mation can be derived directly by considering thermo- 
dynamics in the Matsubara formalism. Eq. (pQ) corre- 
sponds to a trapezoidal-rule approximation of the T = 



integral [23]. At room temperature, the Matsubara fre- 
quency corresponds to a wavelength 27r/£ = 7/im, much 
larger than separations where the Casimir effect is typi- 
cally observed, so usually T > corrections are negligi- 
ble d [l4j . However, experiments are pushing towards 
> 1/im separations [2J-[26| in an attempt to observe 
these corrections. Also, a recent theoretical prediction 
shows much larger T corrections with appropriate mate- 
rial and geometry choices fl6j ] . 

To compute Casimir forces in arbitrary geometries, it 
is desirable to exploit mature methods from computa- 
tional classical electromagnetism (EM), and a number 
of approaches have been suggested [2714311]. One tech- 
nique is to use the fluctuation-dissipation theorem: the 
mean-square electric and magnetic fields (E 2 ) and (H 2 ) 
can thereby be computed from classical Green's functions 
[20| , and the mean stress tensor can be computed and in- 
tegrated to obtain the force p], Q, Q. In particular, at 
each lj, the correlation function of the fields (E 2 ) is given 
by: 



^ Im [co 2 Gf k (u>; x, x')] coth 



(2) 



where G? k = (Gf )j is the classical dyadic "photon" 
Green's function, proportional to the electric field in the 
j direction at x due to an electric-dipole current in the k 
direction at x', and solves 

[V x //(c^x^V x -cj 2 £(cj,x)] Gf (o;,x,x') 

= ^ 3 (x-x , )e fc , (3) 

where e is the electric permittivity tensor, \i is the mag- 
netic permeability tensor, and is a unit vector in direc- 
tion k. The magnetic-field correlation (H 2 ) has a similar 
form p], Q . Note that the temperature dependence ap- 
pears as a coth factor (from a Bose-Einstein distribution). 
If this is Wick-rotated to imaginary frequency lj = i£, the 



2 



poles in the coth function gives the sum (pQ) over Mat- 
subara frequencies [32|. In our EM simulation, what 
is actually computed is the electric or magnetic field in 
response to an electric or magnetic dipole current, re- 
spectively. This is related to Gij by 



(4) 



where Ejk(uo, x, x') denotes the electric field response 
in the jth direction due to a dipole current source 
J(cj, x, x') = 5(x — x')efc 

This equation can be solved for each point on a sur- 
face to integrate the stress tensor, and for each frequency 
to integrate the contributions of fluctuations at all fre- 
quencies. Instead of computing each uo separately, one 
can use a pulse source in time, whose Fourier trans- 
form contains all frequencies. As derived in detail else- 
where p], Q , it turns out that this corresponds to a se- 
quence of time-domain simulations, where pulses of cur- 
rent are injected and some function Y{t) of the resulting 
fields (corresponding to the stress tensor) is integrated 
in time, multiplied by an appropriate weighting factor 
g(t). We perform these simulations by using the stan- 
dard FDTD technique [l?}, which discretizes space and 
time on a uniform grid. In frequency domain, Wick 
rotation to complex cj(£) is crucial for numerical com- 
putations in order to obtain a tractable frequency inte- 
grand d, [23j| , and the analogue in time domain is equally 
important to obtain rapidly decaying fields (and hence 
short simulations) p], [J. In time domain, one must im- 
plement complex uo indirectly: because uj only appears 
explicitly with e in Eq. (j3j), converting uo to the com- 
plex contour oj(£) = £^1 + ^ is equivalent to operat- 
ing at a real frequency _£ with an artificial conductivity 
e(r) —> s(r)(l + y) [l], [2j- One cannot use purely imagi- 
nary frequencies uj = z£ in the time domain, because the 
corresponding material e — » — e has exponentially grow- 
ing solutions in time P. Thus, by adding an artificial 
conductivity everywhere, and including a corresponding 
Jacobian factor in g(t), one obtains the same (physical) 
force result in a much shorter time (with the fields de- 
caying exponentially due to the conductivity). 

Now, we introduce the basic idea of how T > is 
incorporated in the time domain, and explain where the 
difficulty arises. The standard T > analysis of Eq. (pp) 
is expressed in the frequency domain, so we start there 
by exploiting the fact that the time-domain approach is 
derived from a Fourier transform of the frequency-domain 
approach. In particular, g(t) is the Fourier transform 
of a weighting factor in the Fourier domain p], 

0]. At real frequency, the effect of T > is simply to 



include an additional weighting factor coth 
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FIG. 1: Comparison between FDTD (red circles) and the an- 
alytical Lifshitz formula [33|] (blue line) for the Casimir force 
between perfect-metal plates in ID with separation a. The 
uj = and uj ^ contributions to the Matsubara sum JT]) 
are plotted separately, in addition to the total force. The 
straightforward method of including the coth(/kj / '2kT) Boes- 
Einstein factor in the FDTD integration (green dashed line) 
gives an incorrect result because the uj — pole requires spe- 
cial handling. 



naive, approach is to replace with: 



g[u(€)] -+g[u>(Z)]coth 
= -<£ 



ujt 




icr/2£) coth 
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(5) 



integral from Eq. (|2]). So, a straightforward, but 



using the T = g [cj(£)] expression from Ref. Q], and then 
Fourier transform this to yield g(t). However, there is an 
obvious problem with this approach: the 1/uj singularity 

in coth means that Eq. (|5]) is not locally integrable 

around £ = 0, and therefore its Fourier transform is not 
well-defined. If we naively ignore this problem, and com- 
pute the Fourier transform via a discrete Fourier trans- 
form as in [1, 2], simply assigning an arbitrary finite 
value for the £ = term, this unsurprisingly gives an 
incorrect force for T > compared to the analytical Lif- 
shitz formula for the case of parallel perfect-metal plates 
in ID [33|| , as shown in Fig. [T] (green dashed line). 

Instead, a natural solution is to handle uo ^ by the 
coth factor as in Eq. (j5j), but to subtract the uo = pole 
and handle this contribution separately. As explained 
below, we will extract the correct uo = contribution 
from the frequency-domain expression Eq. (pQ), convert it 
to time domain, and add it back in as a manual correction 

to g(t). In particular , the coth function has poles 

at uo = innuoT for integers n. When the frequency integral 
is Wick-rotated to imaginary frequency, the residues of 
these poles give the Matsubara sum Eq. (pQ) via contour 
integration [32|. If we subtract the n = pole from the 
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coth, obtaining 



0n>o(O = 9[w(€)] \ coth 



^(0 



ujt 



(6) 



the result of the time- domain integration of g n> o(t)Y(t) 
will therefore correspond to all of the n > terms in 
Eq. (pp) , nor is there any problem with the Fourier trans- 
formation to g n >o(t)- Precisely this result is shown for 
the ID parallel plates in Fig. [TJ and we see that it indeed 
matches the n > terms from the analytical expression. 
To handle the uj = contribution, we begin with the 
real-u; T = expression for the Casimir force following 
our notation from the time-domain stress-tensor method 



Fj = Im 



duj g R (uj)Ti(uj), 



(7) 



where gR(uj) = — iuo is the weighting factor for the a = 
real uj contour and Yi{uj) = Yf (u) + Yf (uj) is the surface- 
integrated stress tensor (electric- and magnetic-field con- 
tributions). From Eq. (pp), the uj = contribution for 
T > is then 



■ i,(n=0) 



lim Im 

^0+ 



hl ( ■ \V ( ^ 

~o(-^) r iM — * — 

7T Z ft 



lim Re[-ujYi(uj)k B T}. 



(8) 
(9) 



Notice that h cancels in the uj — contribution: 
this term dominates in the limit of large T where the 
fluctuations can be thought of as purely classical ther- 
mal fluctuations. To relate Eq. (|9]) to what is actually 
computed in the FDTD method requires some care be- 
cause of the way in which we transform to the cj(£) 
contour. The quantity Yf(uj) is proportional to an in- 



tegral of Eij(uj) 



jGij(uj), from Eq. (|4|). How- 



ever, the cj(£) transformed system computes ff (£) ~ 
Eij(£) = -^Gij(0, where G(0 solves Eq. © with 
cj 2 e(r) — » £ 2 (1 + y)£(r), but what we actually want is 

— i^Gij^)!^^^) = — ^(0^ij(0- Therefore, the cor- 
rect = contribution is given by 



lim rf (uj) lim 
^0+ £^o+ 



(10) 



Combined with cj(£)/cbT factor from Eq. (j9]), this 
gives an n = contribution of r|^ =0 + multiplied by 
-cj(0 2 &b7Y£|^ + = crk B T. This = term corre- 
sponds to a simple expression in the time domain, since 
r|^ =0 + is simply the time integral of Y(t) and the co- 
efficient ak B T is merely a constant. Therefore, while 
we originally integrated g n> o(t)Y(t) to obtain the n > 
contributions, the n = contribution is included if we 
instead integrate: 



[g n >o(t)+(rk B T]Y(t). 



(11) 



The term [g n >o(t) + o~k B T] generalizes the original g(t) 
function from Ref. 1 to any T > 0. 




FIG. 2: Comparison between FDTD (circles and diamonds) 
and BEM frequency domain (solid and dashed lines) calcu- 
lation of the 2D Casimir force (^-invariant fluctuations) be- 
tween two perfect-metal sidewalls (separation d), normalized 
by the proximate force approximation for the 2D parallel 
plates Fpfa = hc£(3) /8ira 2 . At T = (circles and solid 
lines) total force (black) varies non-monotonically with d, due 
to competition between TE (red) and TM (blue) polarizations 
[2ll ] . At T = 1 x Trch/kBCL (dashed lines and diamonds) BEM 
and FDTD match, but the non-monotonicity disappears. 



We check Eq. ([TT]) for the ID parallel plate case in 
Fig. [1] against the analytical Lifshitz formula [33]. As 
noted above, the g n> o term (|6|) correctly gives the n > 
terms, and we also see that the ak B T term gives the 
correct n = contribution, and hence the total force is 
correct. 

As another check, we consider a more complicated ge- 
ometry: a piston-like configuration from Ref. 4, shown 
schematically in the inset of Fig. [2j This system consists 
of two square rods adjacent between two sidewalls, which 
we solve here for the 2D case of z-invariant fluctuations. 
At T = 0, such geometries were shown to exhibit an in- 
teresting non-monotonic variation of the force between 
the two blocks as a function of sidewall separation d 
0, [H, [2l|, [22I , which does not arise in the simple pairwise- 
interaction heuristic picture of the Casimir force. This 
can be seen in the solid lines of Fig. [2j where the non- 
monotonicity arises from a competition between forces 
from transverse-electric (TE) and transverse-magnetic 
(TM) field polarizations Q, which in turn can be ex- 
plained by a method-of-images argument [2l|. In Fig. [2j 
the solid lines are computed by a T = frequency- 
domain boundary element method (BEM) evaluating a 
path- integral expression |9[ , whereas the circles are com- 
puted by the T = FDTD method [3,[2[, and both meth- 
ods agree. We also compute the force at T = 1 x 7rch/k B a 
where the £ = + term dominates. We see that the 
FDTD method with the T > modification Eq. (JTTJ) (di- 
amonds) agrees with the frequency-domain BEM results 
(dashed lines), where the latter simply use the Matsubara 
sum (PQ) to handle T > 0. 

Interestingly, Fig. [2] shows that the non-monotonic ef- 
fect disappears for T = 1 x i\chjk B a^ despite the fact that 



the method-of-images argument of Ref. HH ostensibly ap- 
plies to the £ = + quasi-static limit (which dominates at 
this large T) as well as to £ > 0. The argument used the 
fact that TM fluctuations can be described by a scalar 
field with Dirichlet boundary conditions (vanishing at the 
metal), and in this case the sidewalls introduce opposite- 
sign mirror sources that reduce the interaction as d de- 
creases; in contrast, TE corresponds to a Neumann scalar 
field (vanishing slope), which requires same-sign mirror 
sources that increase the interaction as d decreases [2l| . 
In Fig. [2j however, while the T = 1 x itch/k^a TM force 
still decreases as d decreases, the TE force no longer in- 
creases for decreasing d at T = 1 x ttcH/ k^,a. The problem 
is that the image-source argument most directly applies 
to z-directed dipole sources in the scalar-field picture — 
electric Jf currents for TM and magnetic currents 
for TE — while the situation for in-plane sources (corre- 
sponding to derivative of the scalar field from dipole-like 
sources) is more complicated [34[ . For a sufficiently large 
T dominated by the £ = + contribution (as is the case 
here), we find numerically that the sources as £ — >• + 
no longer contribute to the force. Intuitively, as £ — >• + 
a magnetic dipole source produces a more and more con- 
stant (long wavelength) field, which automatically satis- 
fies the Neumann boundary conditions and hence is not 
affected by the geometry. Instead, numerical calculations 
show that the TE £ = + contribution is dominated by 
sources and the corresponding electric stress-tensor 
terms, which turn out to slightly decrease in strength as 
d decreases. (A related effect is that, for small it can be 
observed in Fig.[2]that the T = 1 force is actually smaller 
than the T = force, again due to the suppression of the 
TE contribution. Since the force diverges as T — )• oo, this 
means that the force changes non-monotonically with T 
at small d; a similar non-monotonic temperature depen- 



dence was previously observed for Dirichlet scalar-field 
fluctuations in a sphere-plate geometry [l9j.) 

In contrast, if we consider the 3D constant cross- 
section problem with z-dependent fluctuations, corre- 
sponding to integrating e %kzZ fluctuations over k z [20| , 
then we find that the non-monotonic effect is preserved 
at all T. This is easily explained by the fact that, for 
perfect metals, k z ^ is mathe matically equivalent to a 
problem at k z = and £ — >• >/£ 2 + fcf [K, HH, and so 
the n = Matsubara term still contains contributions 
equivalent to £ > in which the mirror argument 
applies and the situation is similar to T = 0. In any 
case, this 2D disappearance of non-monotonicity seems 
unlikely to be experimentally relevant, because we find 
that it only occurs for T > 0.7 x itch/k^a , which for 
a = 1/im corresponds to T > 5000if . 

The main point of this paper is that a simple (but 
not too simple) modification to our previous time-domain 
method allows off-the-shelf FDTD software to easily cal- 
culate Casimir forces at nonzero temperatures. Although 
the disappearance of non-monotonicity employed here as 
a test case appears unrealistic, recent predictions of other 
realistic geometry /material effects [16j, combined with 
the fact that temperature effects in complex geometries 
are almost unexplored at present, lead us to hope that 
future work will reveal further surprising temperature ef- 
fects that are observable in micromechanical systems. 
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